Temperature Relaxation in Hot Dense Hydrogen 
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Temperature equilibration of hydrogen is studied for conditions relevant to inertial confinement 
fusion. New molecular-dynamics simulations and results from quantum many-body theory are 
compared with Landau-Spitzer (LS) predictions for temperatures T from 50 eV to 5000 eV, and 
densities with Wigner-Seitz radii r s = 1.0 and 0.5. The relaxation is slower than the LS result, even 
for temperatures in the keV range, but converges to agreement in the high-T limit. 
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Introduction - While a first-principles description of 
the equilibrium properties of strongly coupled Coulomb 
systems is a formidable task nonequilibrium sys- 
tems pose an even greater challenge. Short-pulse lasers 
and shock waves create nonequilibrium states. Thus, 
Coulomb systems as diverse as warm dense matter 
ultracold plasmas f^, shocked semiconductors and 
dense deuterium [5[ can now be readily created in the 
laboratory, but initially under nonequilibrium conditions. 
Similarly, energy relaxation (ER) in astrophysical plas- 
mas is important to the physics of fusion of H, C, etc., 
in determining stellar evolution^. Here, we consider the 
ER of nonequilibrium dense hydrogen due to its impor- 
tance in inertial confinement fusion (ICF)[7[. 

The earliest theories of ER in plasmas were formulated 
by Landau [|[ and Spitzer (denoted LS). The LS ap- 
proach is applicable to dilute, hot, fully ionized plasms 
where the collisions are weak, binary, and involve neg- 
ligible quantum effects; essentially, LS is Rutherford's 
Coulomb scattering formula applied to a Maxwellian dis- 
tribution. A characteristic feature of the LS approach is 
the use of a Coulomb logarithm (CL), i.e., 



£ = lnA 



db/b- 



dk/k 



(1) 



where b mm (or l/k max ) and b max (or l/k mm ) are suit- 
able, but ad hoc, impact parameter (or momentum) cut- 
offs for the Coulomb collision. The full quantum me- 
chanical method, based on calculating a transition rate, 
does not suffer from this problem. Calculations at the 
Fermi golden-rule level and beyond have been made by 
Dharma-wardana et al. 
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, [HJ . Such methods auto- 
matically include degeneracy effects, effects of collective 
modes, and strong coupling. Other approaches employ 
convergent kinetic equations [H, Hansen and Mc- 
Donald [HI (denoted HM) and Reimann and Toepffer 
have directly obtained ER via molecular dynamics (MD) 
simulations. 

In a hydrogen plasma the particle charges Z\ , Zj 
are ±1, in atomic units, where the electronic charge 
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\e\=h=m e —l. The mean electron and proton densities 
n and p are identical. The ratio of a typical Coulomb 
energy to the kinetic energy becomes, in the classi- 
cal regime, T — l/(r s T), where T is the temperature 
in energy units, and r = r s in the quantum regime. 
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r s = [3/ (47m)] ' is the radius of the Wigner-Seitz sphere 
of an electron or a proton. The properties of partially 
degenerate fully-ionized plasmas require two indepen- 
dent parameters, e.g., both r s and 8 = T/Ep, where 

E F = (3tt 2 7i) 2/3 /2 is the Fermi energy The LS anal- 
ysis of ER is in terms of e-p collisions in a Maxwellian gas. 
Relative to the LS approach, some theoretical approaches 
relax faster [l4[ , or slower [l(| [HJ • HM concluded that 
many-body effects are negligible and that the the LS re- 
sult holds [15j|. There is currently no direct experimental 
data for ER rates, although such experiments are under 
way[H[I!. 

Here, we present larger HM-like molecular-dynamics 
(MD) simulations and quantum many-body calculations 
to narrow the gap in our predictions of ER rates for 
hot, dense hydrogen relevant to ICF targets composed 
of dense cryogenic fuel rapidly laser compressed to keV 
temperatures. The modeling of these dense plasmas cov- 
ers physical processes over many orders of magnitude in 
density and temperature. Nonequilibrium quantum sim- 
ulations require theoretical breakthroughs that are not 
yet fully established; however, the MD techniques that 
employs quantum-corrected effective potentials can be 
usefully applied to these problems. Since MD simula- 
tions attempt to solve the many-body equations of mo- 
tion exactly, they are likely to provide accurate ER rates 
for hot, dense hydrogen. 

Molecular dynamics.- Several issues arise in simula- 
tions of plasmas with temperatures in the 10 2 -10 3 eV 
range. Because the screening length and mean-free- 
path are larger at higher temperatures, we have var- 
ied the number of particles widely (maximum of several 
thousand), with ./V=500 used for the results presented 
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here. Also, because the ER time varies roughly as T e , 
these simulations required millions of time steps for the 
higher-temperature cases. Compounding the longer runs 
was the need for much smaller time steps (as small as 
0.005W" 1 , where oj pe is the usual electron plasma fre- 
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quency), as a result of high velocity collisions at high 
temperature. Such issues are not unexpected, but have 
rarely been dealt with, since MD is typically applied to 
strongly coupled (i.e., cooler) classical systems. We car- 
ried out direct simulations of electrons and protons since 
the Born-Oppcnheimcr approximation is not applicable 
to this problem. The mass ratio (1 : 1836) was used due 
to our interest in mass-dependent collective mode effects. 
Finally, integrations were carried out with the velocity- 
Verlet algorithm using an 0(N 3 / 2 ) Ewald method. 

The Coulomb interaction 1/r of classical physics is re- 
placed by the mean- value of the operator 1/r in quan- 
tum systems. This feature modifies the short-ranged be- 
haviour of the electron-electron and electron-proton in- 
teractions, since the de Broglie wavelength of the elec- 
tron is not negligible for small r. Also, unlike the clas- 
sical electron-proton interaction which always leads to a 
bound state, the delocalized electron does not bind to 
the proton for regimes studied here. These effects are 
included in the e-e and e-p potentials by solving the 
relevant Schrodinger equations for the two-body scat- 
tering processes. These lead to diffraction-corrected 
Coulomb potentials Vy (?"y) which are Coulomb- like for 
distances larger than the respective de Broglie lengths, 
Ay = 1/ ^/27r/iyT e . Here the effective masses of 

the colliding pair. The choice of b m i n = Ay is usal in the 
LS approaches. 

Then interaction </>y(r) between the pair i,j is given 

by 



„ d f r 



(ry) = ~fL 



(2) 

+ (3) 



where 



,Ewl 



is the Ewald potential. The potential, v\ 



accounts for the spin-averaged Pauli exclusion between 
two electrons. This ensures that the "non-interacting" 
electron pair distribution function (PDF) calculated from 
a classical simulation is just the non-interacting quan- 
tum PDF @, EH. The explicit form[22| used by HM 
is adequate for most of the range studied in this paper. 
The factor /y is unity for all except the cross-species 
case f ep . This is chosen using CHNC. The CHNC uses 
the above potentials and a classical fluid temperature 
T c f = \/T^ + Tj where the quantum temperature T q is 
defined in Ref. [2l|. The f ep factor at a given r s , T e , T p is 
fixed by requiring that the CHNC g e p(r), at r = is the 
same as the Kohn-Sham value of g e p(r = 0), as discussed 
more fully in Ref. [23| . The factors f ep are within 5% of 
unity for the conditions of our study. 

To validate the potentials under two-temperature, 
dense-plasma conditions, we carried out (using <^>y and its 
simplified HM forms) MD and CHNC for the conditions 
r s = 1.0, T e = 50eV, and T t = lOeV. In the MD, electron 
and proton velocity-scaling thermostats were used to cre- 
ate the two-temperature system, whereas recent results 
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FIG. 1: Pair distribution functions from MD (data points), 
CHNC (thick lines) for dense hydrogen, r s = 0.5 and 1, 
T/Ef =0.25 and 1. The g ee from a simple HNC calcula- 
tion (i.e., no quantum temperature, no Pauli potential) with 



just vfj r is also shown for r s 



CHNC. The pair distribution functions gy(r) are shown 
in Fig. (HJ). 

Temperature equilibration rates were determined for 
the two densities r s = 0.5 (n e — 1.5 x 10 25 cm -3 ) and r s — 
1.0 (n e = 1.61 x 10 24 cm -3 ) over the temperature range 
0.25 < 9 < 20. In practice, an equilibration stage with 
separate electron and proton velocity-scaling thermostats 
was used to establish a two-temperature system. This 
was followed by a microcanonical evolution in which the 
temperatures were relaxed, with 



N 3 
J i—1 



(4) 



for each species j. Energy conservation was care- 
fully monitored to assure stability at the elevated 
temperatures. Stability was quantified by AE — 

timestep 



Bo 



for the cross-species temperature T e . 



were used in the 



where Ej is the energy at the i-th 
The timestep was chosen by first performing 
several simulations with varying timesteps for temper- 
atures T = 100, 250,500,750eV and noting the impact 
on AE. As mentioned above, the timestep required for 
stability decreases dramatically with increasing tempera- 
ture. A fit to the slope of the temperature profiles yielded 
the equilibration rate; in practice, the ion temperature 
was fixed for all runs at Tj = lOeV so that the reported 



results are dT e /dt. Over the range of temperatures con- 
sidered, the relaxation rate should be very insensitive to 
the ion temperature. 

Quantum transition rates.- The most transparent, 
strictly quantum approach to the calculation of the ER 
rate is to treat it as a transition rate where an electron 
in an initial momentum state ki transfers to a final state 
kf, while a proton in the initial momentum eigenstate 
Pi absorbs energy and transfers to a final eigenstate pf. 
The availability of such states depends on the products 
of Fermi occupation factors 71^(1 — n kf)i an d similarly 
for the proton states. The strength of the transition de- 
pends on the matrix element between the initial and final 
states. This matrix element may be taken in lowest-order 
theory (Born approximation) or in higher order (i.e, a T- 
matrix evaluation). These are the usual ingredients of 
the Fermi golden rule (FGR) for the transition rate. The 
summation over all such pair processes gives the total ER 
rate. But such summations immediately convert the de- 
scription of the plasma into a description containing the 
full spectrum of single-particle and collective modes. The 
spectrum of all modes is given by the spectral function 
Aj(q,ui,Tj) where the species index j = e or p. These 
spectral functions are given by the imaginary parts of the 
corresponding dynamic response functions x J (fc,u;), e.g, 
Eq. (16) of Ref. [II]. The ER rate evaluated within the 
Fermi golden rule, Rf gr can be expressed in terms of the 
response functions of the plasma as follows, given in Eqs. 
(4)-(7) of Ref. [HI, and Eq. (15) of Ref. pj: 



0.2 



TIE. 



R igr 
AB 



dK 
~dt 



d 3 k ujdui 



(2n)" 



2tt 



(AB)F e . 



coth(u/2T e ) - coth(w/2T p ) 



F ep = \(V ep (k)\^ X p (k,u) 



(k,u) 



(5) 

(6) 
(7) 



In the above we have used the spherical symmetry of the 
plasma to write scalars q, k instead of q, k to simplify the 
notation. 

The excess-energy density is denoted by K = K e — Ki, 
and becomes 3(T e — T p )n/2, in the classical regime, i.e., 
where the chemical potential fj. is negative. The interac- 
tion V ep {k) in this equation is the full Coulomb matrix 
element and not the diffraction corrected form used in the 
CHNC and the classical simulations. In the simplest form 
of the FGR, V ep (k) = 4ir/k 2 since the momentum states 
are taken to be plane waves. A T-matrix evaluation 
would use phase-shifted plane waves and the correspond- 
ing modified density of states, instead of d 3 k/ (27r) 3 . 
With the onset of the classical regime where [i < 0, which 
occurs for 9 > 1, the AB factor becomes 2A/u> where 



A = (T e — T p ). It was shown in Hazak et al.[12| how we 
may do the w-integration by exploiting the /-sum rule 
and the fact that the ion-spectral function, peaking near 
the ion-plasma frequency, resides way below the electron 
spectral function. Then Eq. O can be written, to a good 
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FIG. 2: The relaxation time r/1000, in units of the inverse 
electon-plasma frequency, for dense hydrogen, r s = 0.5 and 
1, from the degenerate to the classicle region. The Landau- 
Spitzer (LS) result using the Hansen-McDonald(HM) pre- 
scription for the fc-cutoffs, the results from the Fermi Golden 
rule(FGR), and the coupled-mode(CM) are shown, together 
with the MD results. 



approximation as: 



1 dA 
A~dt 



3 71 



dk (8) 



where uji on is the proton-plasma frequency. If we keep 
the proton temperature T p fixed, we see that Eq. [5] leads 
to a relaxation time r for the electron temperature T e , 
involving the inverse of the r.h.s. of Eq. [5] 

The above analysis treats the plasma as two indepen- 
dent subsystems. In reality, the ion-density fluctuations 
are screened by the electron subsystem, and the ion- 
plasma mode becomes an ion-acoustic mode. The ex- 
citations in the coupled-mode system are described by 
the zeros of Eq. (45) of Ref. [l(|. In the static, k — > 
limit this denominator converts the electron screening 
parameter k e DH to ^J{{k e DH ) 2 + {k p DH ) 2 }. However, the 
proton-density modes act dynamically in the relaxation 
process. Thus the use of static ion screening is incorrect. 
The coupled-mode(CM) approach is fully dynamical and 
includes another denominator, 



4 m = |i-K;(fc)x ee (k,u) x pp (k,u) 



(9) 
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into the integrand in Eq. [5] That is, F ep in Eq. \5\ is 
replaced by F ep /d cm . 

The simple Landau-Spitzer form can also be be writ- 
ten in the same form as Eq. [SJ as shown in Ref. 
The quantum approaches in CM and FGR automatically 
contain the diffraction and screening effects. Thus, while 
Eq. [SJuse the full integration — > oo, LS needs the cut- 
offs k m i n and k max to obtain a convergent result. The 
calculated LS-values of r does depend somewhat on the 
choice of k min and k max . Hence different realizations of 
the LS-form need not reduce to the same result at finite 
T/Ep. In fact Lee and More 24} use cutoffs based on the 



full static screening length which includes the ions as well 
and differ significantly from LS. However, the CM analy- 
sis clearly shows that the ion response in ER is dynamic. 

The non-interacting response function x°(g,w,T) 
at arbitrary degeneracies was given by Khanna and 
Glvde[25| . We use a generalized RPA form where local 
filed corrections G ee (k) may be included [Ttj]. However, 
these are quite small for the conditions of this study. 
Both the FGR and CM calculations assume that linear 
response can be used to discuss the interaction of a pro- 
ton with the electrons. The resulting calculations are 
shown in Fig. [2j The coupled-mode (CM) calculation 
is quite close to the Fermi golden rule (FGR) f-sum re- 
sult. This is expected since the H- plasmas considered 
here are relatively weakly coupled, with V < 1. Never- 



theless, the inclusion of CM leads to better agreement 
with the MD simulation. Also, we have assumed the 
bare 4ir/k 2 form for the V ep in the FGR and CM formu- 
lae, without the moderating effects of a pseudo-potential. 
Such effects would tend to make the t larger than that 
from the present FGR or CM calculation. This linear- 
response assumption is more satisfactory for the r s = 0.5 
plasma. Thus the MD results at r s = 0.5 are very close 
to the CM results. 

Conclusion.- We have evaluated the temperature re- 
laxation time in hot, dense hydrogen using state-of-the- 
art molecular dynamics simulations and quantum many- 
body theory. We find that the relaxation is slower than 
the LS value even for temperatures in the kilovolt range, 
which suggests that burning plasmas are slightly more 
out of equilibrium that might have been expected. Un- 
like in the calculations presented in, say, Ref. [13, EH, 
where strongly coupled Al-plasmas were considered, the 
present calculations are for systems with T ~ 1 or less. 
The temperatures have been pushed to T j Ep ~ 20. Thus 
we see that the CM, FGR and the LS forms converge 
for sufficiently large T/Ep. The MD results are slightly 
higher than from the analytical models which use linear 
response. It is also clear that the LS form is inadequate 
for highly compressed low-T, partially degenerate plas- 
mas. 
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